## Download data
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url = url,
destfile = "GSE63310.tar",
mode = "wb")
utils::untar(tarfile = "GSE63310.tar",
exdir = file.path(tempdir(),"GSE63310"))
## Unzip files
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
"GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
R.utils::gunzip(file.path(tempdir(),"GSE63310/",i), overwrite=TRUE)Análisis de Expresión Diferencial I
Licenciatura en Ciencias Genómicas, Transcriptómica, Semestre 2025-2
1. Motivación
El análisis de expresión diferencial identifica genes o moléculas con cambios significativos en sus niveles de expresión bajo diferentes condiciones, como enfermedades o tratamientos. Es crucial para comprender mecanismos biológicos, descubrir biomarcadores, desarrollar terapias dirigidas y avanzar en la medicina personalizada. En bioinformática, el DEA permite procesar y analizar grandes volúmenes de datos de tecnologías como RNA-Seq o microarrays, ayudando a identificar patrones clave en redes biológicas. Su importancia radica en generar conocimiento sobre enfermedades, guiar investigaciones y acelerar el desarrollo de herramientas diagnósticas y terapéuticas.
1.1 Setup de datos
Vamos a usar este set de datos: https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file
2. EdgeR / Limma
EdgeR es una librería de R que se usa para realizar análisis de expresión diferencial de genes provenientes de datos de RNA-seq. Esta librería emplea métodos estadísticos para experimentos multigrupo basados en modelos lineares generalizados (generalized linear models: glms) los cuales son eficientes para experimentos multifactoriales independientemente de su complejidad. También implementa métodos Bayesianos que permiten la estimación de la variación biológica genética-específica, incluso para experimentos con los mínimos niveles de replicación biológica.
Limma es una librería de R que fue creada para el ánalisis de expresión diferencial para datos provenientes de microarreglos. Usa principalmente modelos lineales para experimentos multigrupos.
2.1 Load libraries
2.2 Read files
## [1] "DGEList"
## attr(,"package")
## [1] "edgeR"
## [1] 27179 9
2.2.1 Reorganize information
## [1] "10_6_5_11" "9_6_5_11" "purep53" "JMS8-2" "JMS8-3" "JMS8-4"
## [7] "JMS8-5" "JMS9-P7c" "JMS9-P8c"
colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
"Basal", "ML", "LP"))
x$samples$group <- group
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
x$samples$lane <- lane
x$samples2.2.2 Add gene names
genes <- AnnotationDbi::select(Mus.musculus,
keys = rownames(x),
columns = c("SYMBOL","TXCHROM"),
keytype = "ENTREZID")
## Remove duplicates
genes <- genes[!duplicated(genes$ENTREZID),]
head(genes)## An object of class "DGEList"
## $samples
## files
## 10_6_5_11 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545535_10_6_5_11.txt
## 9_6_5_11 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545536_9_6_5_11.txt
## purep53 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545538_purep53.txt
## JMS8-2 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545539_JMS8-2.txt
## JMS8-3 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545540_JMS8-3.txt
## JMS8-4 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545541_JMS8-4.txt
## JMS8-5 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545542_JMS8-5.txt
## JMS9-P7c C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545544_JMS9-P7c.txt
## JMS9-P8c C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545545_JMS9-P8c.txt
## group lib.size norm.factors lane
## 10_6_5_11 LP 32863052 1 L004
## 9_6_5_11 ML 35335491 1 L004
## purep53 Basal 57160817 1 L004
## JMS8-2 Basal 51368625 1 L006
## JMS8-3 ML 75795034 1 L006
## JMS8-4 LP 60517657 1 L006
## JMS8-5 Basal 55086324 1 L006
## JMS9-P7c ML 21311068 1 L008
## JMS9-P8c LP 19958838 1 L008
##
## $counts
## Samples
## Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
## 497097 1 2 342 526 3 3 535 2
## 100503874 0 0 5 6 0 0 5 0
## 100038431 0 0 0 0 0 0 1 0
## 19888 0 1 0 0 17 2 0 1
## 20671 1 1 76 40 33 14 98 18
## Samples
## Tags JMS9-P8c
## 497097 0
## 100503874 0
## 100038431 0
## 19888 0
## 20671 8
## 27174 more rows ...
##
## $genes
## ENTREZID SYMBOL TXCHROM
## 1 497097 Xkr4 chr1
## 2 100503874 Gm19938 <NA>
## 3 100038431 Gm10568 <NA>
## 4 19888 Rp1 chr1
## 5 20671 Sox17 chr1
## 27174 more rows ...
2.3 Data preprocessing
## counts per million
cpm.obj <- cpm(x)
## log counts per million transformation
lcpm.obj <- cpm(x, log = TRUE)
## Mean and median
L <- mean(x$samples$lib.size) * 1e-6
M <- median(x$samples$lib.size) * 1e-6
c(L, M)## [1] 45.48855 51.36862
## The average library size is about 45.5 million and the minimum
## log-CPM should be log2(2/45.5) = -4.51 which means that
## a count of 0 maps to a log-CPM value of -4.51
summary(lcpm.obj)## 10_6_5_11 9_6_5_11 purep53 JMS8-2
## Min. :-4.5074 Min. :-4.5074 Min. :-4.50743 Min. :-4.5074
## 1st Qu.:-4.5074 1st Qu.:-4.5074 1st Qu.:-4.50743 1st Qu.:-4.5074
## Median :-0.6847 Median :-0.3589 Median :-0.09513 Median :-0.0901
## Mean : 0.1714 Mean : 0.3312 Mean : 0.43559 Mean : 0.4089
## 3rd Qu.: 4.2913 3rd Qu.: 4.5601 3rd Qu.: 4.60081 3rd Qu.: 4.5475
## Max. :14.7632 Max. :13.4952 Max. :12.95700 Max. :12.8513
## JMS8-3 JMS8-4 JMS8-5 JMS9-P7c
## Min. :-4.5074 Min. :-4.5074 Min. :-4.50743 Min. :-4.5074
## 1st Qu.:-4.5074 1st Qu.:-4.5074 1st Qu.:-4.50743 1st Qu.:-4.5074
## Median :-0.4281 Median :-0.4064 Median :-0.07152 Median :-0.1704
## Mean : 0.3225 Mean : 0.2529 Mean : 0.40428 Mean : 0.3708
## 3rd Qu.: 4.5772 3rd Qu.: 4.3199 3rd Qu.: 4.42513 3rd Qu.: 4.6031
## Max. :12.9578 Max. :14.8520 Max. :13.19491 Max. :12.9413
## JMS9-P8c
## Min. :-4.5074
## 1st Qu.:-4.5074
## Median :-0.3300
## Mean : 0.2749
## 3rd Qu.: 4.4355
## Max. :14.0102
2.4 Remove genes that are lowly expressed
2.5 Normalize gene expression
Se necesita normalizar debido a que los datos de RNA-Seq contienen sesgos técnicos y biológicos que pueden distorsionar los resultados si no se corrigen. Las razones principales para normalizar son: diferencias de profundidad de secuenciación, diferencias en la composición génica, comparabilidad entre muestras, impacto de la longitud de los genes, preparación para métodos estadísticos.
Métodos comunes de normalización:
- TMM (Trimmed Mean of M-values): Corrige la profundidad de secuenciación y la composición génica.
- RPKM/FPKM (Reads/Fragments Per Kilobase Million): Normaliza por la longitud del gen y el tamaño total de las muestras.
- TPM (Transcripts Per Million): Similar a RPKM, pero mejora la comparabilidad entre muestras.
- DESeq2/EdgeR Scaling Factors: Ajustan los conteos para análisis estadísticos robustos.
## Normalize
x <- calcNormFactors(x,method = "TMM") # Trimmed mean of M-values
x$samples$norm.factors## [1] 0.8943956 1.0250186 1.0459005 1.0458455 1.0162707 0.9217132 0.9961959
## [8] 1.0861026 0.9839203
Comparación entre normalizar o no
Visualización:
col <- c("#A6CEE3","#1F78B4","#B2DF8A","#33A02C","#FB9A99","#E31A1C","#FDBF6F","#FF7F00","#CAB2D6")
par(mfrow=c(1,2))
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="A. Example: Unnormalised data",ylab="Log-cpm")
x2 <- calcNormFactors(x2)
x2$samples$norm.factors## [1] 0.05770899 6.08287835 1.22023972 1.16478991 1.19661094 1.04659233 1.15048074
## [8] 1.25431164 1.10901983
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="B. Example: Normalised data",ylab="Log-cpm")2.6 Unsupervised clustering of samples
- Usar una MultiDimensional Scaling plot (MDS): esta muestra la similitudes y disimilitudes entre las muestras con algoritmo no supervisado. La primera dimensión representa el fold-change que mejor separa a las muestras y que explica la mayor proporción de la variancia de los datos. Aquí nos podemos dar cuenta si se deben de hacer algunas correcciones del batch.
lcpm <- cpm(x, log=TRUE)
## Define colors
par(mfrow=c(1,2))
col.group <- group
levels(col.group) <- RColorBrewer::brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)
## Define colors
col.lane <- lane
levels(col.lane) <- RColorBrewer::brewer.pal(nlevels(col.lane), "Set2")
col.lane <- as.character(col.lane)
## Plot MDS
plotMDS(lcpm, labels=group, col=col.group)
title(main="A. Sample groups")
plotMDS(lcpm, labels=lane, col=col.lane, dim=c(3,4))
title(main="B. Sequencing lanes")2.7 Differential expression analysis
2.7.1 Create a matrix
Genes que están expresados en diferentes niveles. Los modelos lineales se ajustan a los datos bajo la asunción de que los datos siguen una distribución normal (lo cual no es cierto porque siguen una distribución binomial negativa o una de poisson).
Este diseño experimental quita el intercepto del primer factor group pero un intercepto permanece con el segundo factor lane. Es esencial entender cómo interpretar los coeficientes en el modelo.
design <- model.matrix(~0 + group + lane)
colnames(design) <- gsub("group","",colnames(design))
design## Basal LP ML laneL006 laneL008
## 1 0 1 0 0 0
## 2 0 0 1 0 0
## 3 1 0 0 0 0
## 4 1 0 0 1 0
## 5 0 0 1 1 0
## 6 0 1 0 1 0
## 7 1 0 0 1 0
## 8 0 0 1 0 1
## 9 0 1 0 0 1
## attr(,"assign")
## [1] 1 1 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
##
## attr(,"contrasts")$lane
## [1] "contr.treatment"
contr.matrix <- makeContrasts(
BasalvsLP = Basal-LP,
BasalvsML = Basal - ML,
LPvsML = LP - ML,
levels = colnames(design))
contr.matrix## Contrasts
## Levels BasalvsLP BasalvsML LPvsML
## Basal 1 1 0
## LP -1 0 1
## ML 0 -1 -1
## laneL006 0 0 0
## laneL008 0 0 0
2.7.2 Removing heteroscedascity from count data
Para los datos de RNA-seq la varianza no es independiente de la media sin importar si los datos sufren una transformación logarítmica o no. Los modelos que asumen una distribución binomial negativa asumen una relación cuadrática media-varianza.
El plot de voom ilustra la relación media-varianza. Este plot muestra una trend que decrece entre la media y la varianza como resultado de la combinación de la varianza técnica en los experimentos de secuenciación y de la varianza biológica entre las réplicas de las muestras de las diferentes poblaciones celulares.
Experimentos con una variación biológica alta resultan en tendencias más planas.
Experimentos con una variación biológica baja reusltan en tendencias decrecientes pronunciadas.
Si no se tienen más de 3 réplicas se opta por usar un GLM o si se desea usar un modelo específico para RNA-seq.
## An object of class "EList"
## $genes
## ENTREZID SYMBOL TXCHROM
## 1 497097 Xkr4 chr1
## 5 20671 Sox17 chr1
## 6 27395 Mrpl15 chr1
## 7 18777 Lypla1 chr1
## 9 21399 Tcea1 chr1
## 16619 more rows ...
##
## $targets
## files
## 10_6_5_11 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545535_10_6_5_11.txt
## 9_6_5_11 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545536_9_6_5_11.txt
## purep53 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545538_purep53.txt
## JMS8-2 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545539_JMS8-2.txt
## JMS8-3 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545540_JMS8-3.txt
## JMS8-4 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545541_JMS8-4.txt
## JMS8-5 C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545542_JMS8-5.txt
## JMS9-P7c C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545544_JMS9-P7c.txt
## JMS9-P8c C:\\Users\\chomb\\AppData\\Local\\Temp\\RtmpwZFj5a/GSE63310//GSM1545545_JMS9-P8c.txt
## group lib.size norm.factors lane
## 10_6_5_11 LP 29387429 0.8943956 L004
## 9_6_5_11 ML 36212498 1.0250186 L004
## purep53 Basal 59771061 1.0459005 L004
## JMS8-2 Basal 53711278 1.0458455 L006
## JMS8-3 ML 77015912 1.0162707 L006
## JMS8-4 LP 55769890 0.9217132 L006
## JMS8-5 Basal 54863512 0.9961959 L006
## JMS9-P7c ML 23139691 1.0861026 L008
## JMS9-P8c LP 19634459 0.9839203 L008
##
## $E
## Samples
## Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5
## 497097 -4.292165 -3.856488 2.5185849 3.2931366 -4.459730 -3.994060 3.2869677
## 20671 -4.292165 -4.593453 0.3560126 -0.4073032 -1.200995 -1.943434 0.8442767
## 27395 3.876089 4.413107 4.5170045 4.5617546 4.344401 3.786363 3.8990635
## 18777 4.708774 5.571872 5.3964008 5.1623650 5.649355 5.081611 5.0602470
## 21399 4.785541 4.754537 5.3703795 5.1220551 4.869586 4.943840 5.1384776
## Samples
## Tags JMS9-P7c JMS9-P8c
## 497097 -3.2103696 -5.295316
## 20671 -0.3228444 -1.207853
## 27395 4.3396075 4.124644
## 18777 5.7513694 5.142436
## 21399 5.0308985 4.979644
## 16619 more rows ...
##
## $weights
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1.079413 1.332986 19.826915 20.273253 1.993686 1.395853 20.494977
## [2,] 1.170357 1.456380 4.804866 8.660025 3.612508 2.626870 8.760149
## [3,] 20.219073 25.573792 30.434759 28.528310 31.352260 25.743247 28.722497
## [4,] 26.947557 32.505933 33.583128 33.232125 34.231754 32.354158 33.334340
## [5,] 26.610864 28.501638 33.645479 33.206374 33.573492 31.996623 33.308490
## [,8] [,9]
## [1,] 1.107780 1.079413
## [2,] 3.211473 2.541942
## [3,] 21.200072 16.657930
## [4,] 30.348630 24.259801
## [5,] 25.171513 23.573305
## 16619 more rows ...
##
## $design
## Basal LP ML laneL006 laneL008
## 1 0 1 0 0 0
## 2 0 0 1 0 0
## 3 1 0 0 0 0
## 4 1 0 0 1 0
## 5 0 0 1 1 0
## 6 0 1 0 1 0
## 7 1 0 0 1 0
## 8 0 0 1 0 1
## 9 0 1 0 0 1
## attr(,"assign")
## [1] 1 1 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
##
## attr(,"contrasts")$lane
## [1] "contr.treatment"
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")Si no se tienen más de 3 réplicas se opta por usar un GLM o si se desea usar un modelo específico para RNA-seq.
## Adjust GLM
x <- estimateDisp(x, design)
fit <- glmQLFit(x, design)
## Make contrast matrices
contr.matrix <- makeContrasts(
BasalvsLP = Basal-LP,
BasalvsML = Basal - ML,
LPvsML = LP - ML,
levels = colnames(design))
contr.matrix## Contrasts
## Levels BasalvsLP BasalvsML LPvsML
## Basal 1 1 0
## LP -1 0 1
## ML 0 -1 -1
## laneL006 0 0 0
## laneL008 0 0 0
## Fit the model
qlf_basal_vs_lp <- glmQLFTest(fit, contrast=contr.matrix[,"BasalvsLP"])
qlf_basal_vs_ml <- glmQLFTest(fit, contrast=contr.matrix[,"BasalvsML"])
## Results
top_basal_vs_lp <- topTags(qlf_basal_vs_lp, n=Inf)
head(top_basal_vs_lp$table)## 1*Basal -1*LP
## Down 4613
## NotSig 7032
## Up 4979
## 1*Basal -1*ML
## Down 4889
## NotSig 6852
## Up 4883
2.7.1 Examining the number of DE genes
La significancia está dada por el threshold al p-value ajustado (5%). Aunque a veces se emplea el valor del log-foldchange como threshold.
## BasalvsLP BasalvsML LPvsML
## Down 4648 4927 3135
## NotSig 7113 7026 10972
## Up 4863 4671 2517
## BasalvsLP BasalvsML LPvsML
## Down 1632 1777 224
## NotSig 12976 12790 16210
## Up 2016 2057 190
## Extract the DEG's shared between both conditions
de.common <- which(dt[,1]!=0 & dt[,2]!=0)
length(de.common)## [1] 2784
2.7.1 Examining individual DE genes from top to bottom
Extraer los resultados más significativos.
2.7.1 Useful graphical representations of differential expression results
MD plot
Volcano plot
par(mfrow=c(1,2))
plotVolcano(df = basal.vs.lp,
pval_col = "P.Value",
lfc_col = "logFC",
adjpval_col = "adj.P.Val",
plot_title = "Basal vs. LP")## [1] TRUE
plotVolcano(df = basal.vs.ml,
pval_col = "P.Value",
lfc_col = "logFC",
adjpval_col = "adj.P.Val",
plot_title = "Basal vs. LP")## [1] TRUE
Heatmap
## Basal vs. LP heatmap
basal.vs.lp.topgenes <- basal.vs.lp$ENTREZID[1:100]
i <- which(v$genes$ENTREZID %in% basal.vs.lp.topgenes)
hmap.palette <- colorRampPalette(c("red", "white", "blue"))
pheatmap::pheatmap(lcpm[i,],
scale = "row",
cluster_rows = TRUE,
labels_row = v$genes$SYMBOL[i],
labels_col = group,
color = hmap.palette(100),
main = "Basal vs. LP")## Basal vs. ML heatmap
basal.vs.ml.topgenes <- basal.vs.ml$ENTREZID[1:100]
i <- which(v$genes$ENTREZID %in% basal.vs.ml.topgenes)
hmap.palette <- colorRampPalette(c("red", "white", "blue"))
pheatmap::pheatmap(lcpm[i,],
scale = "row",
cluster_rows = TRUE,
labels_row = v$genes$SYMBOL[i],
labels_col = group,
color = hmap.palette(100),
main = "Basal vs. ML")3. DESeq2
DESeq2 es ubna librería de R que sirve para realizar experimentos de expresión diferencial. Su input es la matriz de conteos crudos resultantes de un experimento de secuenciación de RNA-seq.
3.1 Load libraries
3.2 Read files
files <- file.path(tempdir(),"GSE63310/",files)
## Load all the files into a list of dataframes
counts.list <- lapply(files, function(file) {
data <- read.delim(file, header = TRUE, stringsAsFactors = FALSE)
data[, c("EntrezID", "Count")]
})
names(counts.list) <- substring(basename(files), 12, nchar(basename(files)) - 4)
## Merge the dataframes into a single counts matrix
counts.mtx <- Reduce(function(x, y) {
merge(x, y, by = "EntrezID", all = TRUE)
}, counts.list)
rownames(counts.mtx) <- counts.mtx$EntrezID
counts.mtx <- counts.mtx[, -1]
colnames(counts.mtx) <- names(counts.list)3.2.1 Create metadata
Es neceario que se cree un dataframe con los metadatos del experimento. Aquí una muestra de las condiciones experimentales.
Group: condiciones experimentales
Lane: se refiere al batch (lote, ronda de secuenciación)
lib.size: se trata del total de counts registrados
norm.factors: factores de normalización
3.3 Create DESeq2 object
Se puede importar de diferentes fuentes. Una de ellas es usando archivos de la abundancia de transcritos con la librería tximport, otra es con una matriz de conteos, otra con archivos de conteos de htseq-count y la última de un objeto de SummarizedExperiment.
## Create DESeq2 object
dds <- DESeqDataSetFromMatrix(countData = counts.mtx,
colData = dds_metadata,
design = ~ group + lane)
dds## class: DESeqDataSet
## dim: 27179 9
## metadata(1): version
## assays(1): counts
## rownames(27179): 11287 11298 ... 100862405 100862406
## rowData names(0):
## colnames(9): 10_6_5_11 9_6_5_11 ... JMS9-P7c JMS9-P8c
## colData names(4): group lib.size norm.factors lane
3.4 Data pre-filter
Se necesita realizar un pre filtrado de los genes que tuvieron un conteo bajo. Por un lado se reduce el tamaño que ocupa en la memoria el objeto dds e incrementar la velocidad de modelado de los conteos. Como valores standard se ocupa lo siguiente:
rows con un mínimo de 10 counts.
Usar el smallest group size el cual es igual al número de muestras.
## We have 9 samples
smallestGroupSize <- 9
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep,]
dds## class: DESeqDataSet
## dim: 13274 9
## metadata(1): version
## assays(1): counts
## rownames(13274): 11302 11303 ... 100862307 100862400
## rowData names(0):
## colnames(9): 10_6_5_11 9_6_5_11 ... JMS9-P7c JMS9-P8c
## colData names(4): group lib.size norm.factors lane
3.4.1 Refactor metadata
Adicionalmente se pueden colapsar todas las réplicas en columnas individuales con la función DESeq2::collapseReplicates() ESTO SOLO PARA RÉPLICAS TÉCNICAS, NO PARA BIOLÓGICAS (una réplica técnica implica múltiples corridas de secuenciación de la misma librería).
3.4.2 Add gene names
3.5 Differential Expression Analysis
Así de simple.
## log2 fold change (MLE): lane L008 vs L004
## Wald test p-value: lane L008 vs L004
## DataFrame with 13274 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## 11302 524.8568 0.2025970 0.427536 0.473871 0.635592 0.914339
## 11303 5688.7471 0.1132199 0.312441 0.362372 0.717074 0.935875
## 11304 43.6222 0.1260263 0.510059 0.247082 0.804845 0.960693
## 11305 3258.0083 -0.0949987 0.168277 -0.564539 0.572387 0.893517
## 11306 1850.6776 0.0372121 0.181120 0.205455 0.837216 0.967711
## ... ... ... ... ... ... ...
## 100862251 1084.2679 0.157915 0.199627 0.791051 0.428914 0.842320
## 100862257 3711.2564 -0.180482 0.157793 -1.143787 0.252712 0.729788
## 100862258 67.0328 -0.518082 0.381632 -1.357545 0.174608 0.646332
## 100862307 51.7442 0.321677 0.394997 0.814379 0.415428 0.834233
## 100862400 88.7577 0.359888 0.497736 0.723049 0.469650 0.857534
Si quisiéramos extraer los genes diferencialmente expresados filtrando por el adjusted p-value entonces se le debe de agregar el parámetro alpha = <float> a la función results(). Por ejemplo, si queremos los DEG’s que son menores o iguales a 0.05 entonces se ve así:
##
## out of 13274 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 112, 0.84%
## LFC < 0 (down) : 196, 1.5%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 17)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## [1] 308
Además, se pueden revisar las comparaciones hechas:
## [1] "Intercept" "group_LP_vs_Basal" "group_ML_vs_Basal"
## [4] "lane_L006_vs_L004" "lane_L008_vs_L004"
3.5.1 Exploring results
Para la visualización se pueden “encoger” los resultados de acuerdo al efecto del tamaño. El estimador por default es apeglm el cual estima con una tendencia al 0, lo que promueve la estabilidad de los estimados en los cuales hay una variabilidad alta o tamaños de muestra pequeños.
MA-plot
Get dispersion factors
Aquí se obtienen los factores de dispersión que es similar a lo que hace la función voom en EdgeR. Aquí se usa la función para Variance stabilizing transformations (vst) la cual transforma a los datos en una escala log2 la cual ha sido normalizada con respecto al tamaño de la librería u otros factores de normalización. La finalidad de esto es eliminar la dependencia de la varianza en la media (si se trabaja con niveles de expresión, entonces genes con baja expresión tendrán una desviación estándar mayor que los genes de alta expresión debido al tipo de distribución de los datos de RNA-seq). El parámetro blind es TRUE por default porque este re-estima las dispersiones empleando solo un intercepto. Esto no es adecuado si se espera que la mayoría de los genes tengan grandes diferencias en sus conteos (lo que puede causar grandes estimados de dispersión por el diseño experimental y no por una diferencia biológica).
3.5.2 Extract the top DEG’s from top to bottom
Una vez que ya se tiene la normalización con vsd, se puede proceder a extraer los DEG’s de interés para cada una de las comparaciones.
## Extract the results for each condition
res_LP_vs_Basal <- results(dds, contrast = c("group", "LP", "Basal"))
res_ML_vs_Basal <- results(dds, contrast = c("group", "ML", "Basal"))
## Add gene names
res_LP_vs_Basal$SYMBOL <- rowData(dds[rownames(res_LP_vs_Basal)])$SYMBOL
res_ML_vs_Basal$SYMBOL <- rowData(dds[rownames(res_ML_vs_Basal)])$SYMBOL
res_LP_vs_Basal$ENTREZID <- rowData(dds[rownames(res_LP_vs_Basal)])$ENTREZID
res_ML_vs_Basal$ENTREZID <- rowData(dds[rownames(res_ML_vs_Basal)])$ENTREZID
## Filter the results
res_LP_vs_Basal_sig <- res_LP_vs_Basal[which(res_LP_vs_Basal$padj < 0.05), ]
res_LP_vs_Basal_sig <- res_LP_vs_Basal_sig[order(res_LP_vs_Basal_sig$padj, decreasing = FALSE), ]
res_ML_vs_Basal_sig <- res_ML_vs_Basal[which(res_ML_vs_Basal$padj < 0.05), ]
res_ML_vs_Basal_sig <- res_ML_vs_Basal_sig[order(res_ML_vs_Basal_sig$padj, decreasing = FALSE), ]Volcano plot
par(mfrow=c(1,2))
plotVolcano(df = res_LP_vs_Basal_sig,
pval_col = "pvalue",
lfc_col = "log2FoldChange",
adjpval_col = "padj",
plot_title = "Basal vs. LP")## [1] TRUE
plotVolcano(df = res_ML_vs_Basal_sig,
pval_col = "pvalue",
lfc_col = "log2FoldChange",
adjpval_col = "padj",
plot_title = "Basal vs. ML")## [1] TRUE
Heatmap
## Basal vs. LP heatmap
basal.vs.lp.topgenes.deseq <- res_LP_vs_Basal_sig$ENTREZID[1:100]
i <- which(rowData(dds)$ENTREZID %in% basal.vs.lp.topgenes.deseq)
hmap.palette <- colorRampPalette(c("red", "white", "blue"))
pheatmap::pheatmap(assay(vsd)[i,],
scale = "row",
cluster_rows = TRUE,
labels_row = rowData(dds[i])$SYMBOL,
labels_col = group,
color = hmap.palette(100),
main = "Basal vs. LP")## Basal vs. ML heatmap
basal.vs.ml.topgenes.deseq <- res_ML_vs_Basal_sig$ENTREZID[1:100]
i <- which(rowData(dds)$ENTREZID %in% basal.vs.ml.topgenes.deseq)
hmap.palette <- colorRampPalette(c("red", "white", "blue"))
pheatmap::pheatmap(assay(vsd)[i,],
scale = "row",
cluster_rows = TRUE,
labels_row = rowData(dds[i])$SYMBOL,
labels_col = group,
color = hmap.palette(100),
main = "Basal vs. ML")4. pyDESeq2
Este es una librería de Pyhton la cual pretende implementar la librería original de R DESeq2 para el análisis de expresión diferencial con datos de bulk RNA-seq. Como es una implementación desde 0, hay algunos cambios importantes en términos de las funciones y de los valores que estas regresan.
4.1 Load libraries
import re
import numpy as np
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
# %pip install sanbomics
from sanbomics.plots import volcano
from sanbomics.tools import id_map
# %pip install pydeseq2
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference
from pydeseq2.ds import DeseqStats
# %pip install scanpy
import scanpy as sc4.2 Read files
## Import R files variable to python
files = r.files
## Empty df
dfs = []
## Read files
for file in files:
df = pd.read_csv(file, sep="\t",usecols = ["EntrezID","Count"])
column_name = file.split("/")[-1].replace(".txt","")
column_name = re.sub(r"^GSM15455\d{2}_", "", column_name)
df.rename(columns={"Count": column_name}, inplace=True)
dfs.append(df)
## Create the matrix
merged_matrix = pd.concat(
[df.set_index("EntrezID") for df in dfs], axis=1
).fillna(0)4.2.1 Create metadata
4.2.2 Add genes
4.2.3 Pre-filter counts
4.3 Create pyDESeq2 object
4.4 Differential Expression Analysis
## Fitting size factors...
## ... done in 0.00 seconds.
##
## Fitting dispersions...
## ... done in 1.28 seconds.
##
## Fitting dispersion trend curve...
## ... done in 0.15 seconds.
##
## Fitting MAP dispersions...
## ... done in 1.35 seconds.
##
## Fitting LFCs...
## ... done in 0.69 seconds.
##
## Calculating cook's distance...
## ... done in 0.00 seconds.
##
## Replacing 0 outlier genes.
## Retrieve stats
result = DeseqStats(dds,inference = inference)
basal_vs_lp = DeseqStats(dds, contrast=["group", "Basal", "LP"], inference=inference)
basal_vs_ml = DeseqStats(dds, contrast=["group", "Basal", "ML"], inference=inference)
## Extract the results
result.summary()## Running Wald tests...
## ... done in 0.38 seconds.
##
## Log2 fold change & Wald test p-value: lane L006 vs L004
## baseMean log2FoldChange lfcSE stat pvalue padj
## 11302 782.020723 -0.123189 0.202669 -0.607834 0.543297 0.875635
## 11303 1651.935164 0.039150 0.159773 0.245032 0.806431 0.959511
## 11304 1344.515458 0.006769 0.135751 0.049865 0.960230 0.992764
## 11305 295.742823 -0.399928 0.211048 -1.894967 0.058097 0.416339
## 11306 2356.502410 0.055302 0.137805 0.401304 0.688196 0.926800
## ... ... ... ... ... ... ...
## 100862251 838.908288 0.297839 0.232688 1.279994 0.200547 0.660890
## 100862257 1041.025804 -0.139267 0.196478 -0.708816 0.478438 0.843321
## 100862258 1534.037447 -0.227660 0.191033 -1.191730 0.233367 0.692072
## 100862307 369.536157 -0.355791 0.216680 -1.642006 0.100589 0.523408
## 100862400 855.267041 -0.344427 0.208666 -1.650609 0.098818 0.520729
##
## [13274 rows x 6 columns]
## Running Wald tests...
## ... done in 0.34 seconds.
##
## Log2 fold change & Wald test p-value: group Basal vs LP
## baseMean log2FoldChange ... pvalue padj
## 11302 782.020723 0.523221 ... 2.276141e-02 4.015082e-02
## 11303 1651.935164 0.319091 ... 7.795831e-02 1.195217e-01
## 11304 1344.515458 0.397116 ... 9.764703e-03 1.865793e-02
## 11305 295.742823 -5.521131 ... 1.947981e-77 1.213967e-75
## 11306 2356.502410 -0.791075 ... 3.972588e-07 1.463156e-06
## ... ... ... ... ... ...
## 100862251 838.908288 0.128801 ... 6.241817e-01 6.960756e-01
## 100862257 1041.025804 0.187428 ... 4.008521e-01 4.868171e-01
## 100862258 1534.037447 0.143802 ... 5.059068e-01 5.888641e-01
## 100862307 369.536157 0.047139 ... 8.487170e-01 8.857512e-01
## 100862400 855.267041 -0.072764 ... 7.580740e-01 8.120299e-01
##
## [13274 rows x 6 columns]
## Running Wald tests...
## ... done in 0.35 seconds.
##
## Log2 fold change & Wald test p-value: group Basal vs ML
## baseMean log2FoldChange ... pvalue padj
## 11302 782.020723 0.125763 ... 5.822965e-01 6.564249e-01
## 11303 1651.935164 -0.310793 ... 8.461658e-02 1.275209e-01
## 11304 1344.515458 0.460100 ... 2.679736e-03 5.633643e-03
## 11305 295.742823 -5.023438 ... 2.219609e-64 8.299461e-63
## 11306 2356.502410 -0.979360 ... 3.187953e-10 1.490556e-09
## ... ... ... ... ... ...
## 100862251 838.908288 0.708760 ... 7.059043e-03 1.372517e-02
## 100862257 1041.025804 -1.106419 ... 5.835175e-07 2.015512e-06
## 100862258 1534.037447 0.310210 ... 1.509529e-01 2.109996e-01
## 100862307 369.536157 -0.459937 ... 6.060102e-02 9.521993e-02
## 100862400 855.267041 0.166803 ... 4.797587e-01 5.609369e-01
##
## [13274 rows x 6 columns]
4.4.1 Add gene names
## Add the genes information
res = res.merge(dds.var,left_index=True,right_index=True)
basal_vs_lp_res = basal_vs_lp_res.merge(dds.var,left_index=True,right_index=True)
basal_vs_ml_res = basal_vs_ml_res.merge(dds.var,left_index=True,right_index=True)
## Sort by adjusted p-value
basal_vs_lp_res = basal_vs_lp_res.sort_values(by='padj', ascending=True)
basal_vs_ml_res = basal_vs_ml_res.sort_values(by='padj', ascending=True)4.4.2 Exploring results
MA-plot
## MA-plot
plt.figure(figsize=(8, 6))
sns.scatterplot(
x=res["baseMean"],
y=res["log2FoldChange"],
hue=res["padj"] < 0.05, # Highlight significant points
palette={True: "blue", False: "gray"},
alpha=0.6,
)
plt.axhline(0, color="black", linestyle="--", linewidth=1)
plt.xscale("log")
plt.xlabel("Mean Normalized Counts")
plt.ylabel("Log2 Fold Change")
plt.title("MA Plot")
plt.legend(title="Significant", labels=["Yes", "No"])
plt.tight_layout()
plt.show()PCA
sc.tl.pca(dds)
## Group PCA
sc.pl.pca(dds, color='group', size=200, annotate_var_explained=True,legend_loc='best', show=True)## Lane PCA
sc.pl.pca(dds, color='lane', size=200, annotate_var_explained=True,legend_loc='best', show=True)Volcano plot
## Basal vs. LP volcano plot
volcano(
basal_vs_lp_res,
symbol = 'SYMBOL',
to_label= 30,
pval_thresh = 0.05,
log2fc_thresh = 2.5)## 0s encountered for p value, imputing 1e-323
## impute your own value if you want to avoid this
## Basal vs. LP volcano plot
volcano(
basal_vs_lp_res,
symbol = 'SYMBOL',
to_label= 30,
pval_thresh = 0.05,
log2fc_thresh = 2.5)## 0s encountered for p value, imputing 1e-323
## impute your own value if you want to avoid this
Heatmap
## Extract the top genes
basal_vs_lp_topgenes = basal_vs_lp_res.ENTREZID[0:100]
basal_vs_ml_topgenes = basal_vs_ml_res.ENTREZID[0:100]
## Add log transformation counts
dds.layers['log1p'] = np.log1p(dds.layers['normed_counts'])
## Get index from significant genes
basalvslp_sigs = dds[:, basal_vs_lp_topgenes]
basalvsml_sigs = dds[:, basal_vs_ml_topgenes]
## Create counts matrices
grapher_blp = pd.DataFrame(
basalvslp_sigs.layers['log1p'].T,
index=basalvslp_sigs.var_names,
columns=basalvslp_sigs.obs_names)
grapher_bml = pd.DataFrame(
basalvsml_sigs.layers['log1p'].T,
index=basalvsml_sigs.var_names,
columns=basalvsml_sigs.obs_names)
## Define color palette
hmap_palette = sns.diverging_palette(10, 250, as_cmap=True)
## Heatmap for Basal vs. LP
sns.clustermap(grapher_blp,
z_score=0,
cmap=hmap_palette,
col_cluster=True,
row_cluster=True,
xticklabels=metadata.group[grapher_bml.columns].to_list(),
yticklabels=basal_vs_lp_res.SYMBOL[basal_vs_lp_topgenes].to_list()).fig.suptitle('Basal vs. LP')
plt.show()## Heatmap for Basal vs. ML
sns.clustermap(grapher_bml,
z_score=0,
cmap=hmap_palette,
col_cluster=True,
row_cluster=True,
xticklabels=metadata.group[grapher_bml.columns],
yticklabels=basal_vs_ml_res.SYMBOL[basal_vs_ml_topgenes].to_list()).fig.suptitle('Basal vs. ML')
plt.show()5. Comparación
5.1 edgeR
Método principal:
Utiliza modelos lineales generalizados (GLM) basados en la distribución binomial negativa.Normalización:
Realiza la normalización utilizando el método de la media recortada ponderada de la razón (TMM). Este método ajusta los efectos de composición entre muestras.Estadísticas diferenciales:
Estima la dispersión biológica para cada gen y utiliza un enfoque Bayesiano empírico para estabilizar estas estimaciones, mejorando la precisión en experimentos con pocos replicados.Casos de uso:
Ideal para datos con pocas réplicas biológicas.
Muy robusto frente a variaciones en la profundidad de secuenciación.
5.2 DESeq2
Método principal:
Basado en la distribución binomial negativa, modela los datos de conteo utilizando un enfoque de regresión GLM, similar a edgeR. Sin embargo, DESeq2 implementa un procedimiento más detallado para la estimación de dispersión.Normalización:
Utiliza el método de razón geométrica para calcular factores de escala, que corrigen las diferencias de tamaño de biblioteca entre muestras.Estadísticas diferenciales:
Estima la dispersión utilizando un modelo empírico y luego ajusta estos valores utilizando un enfoque Bayesiano para estabilizar las varianzas.
Implementa un procedimiento de “shrinkage” para los coeficientes del modelo, lo que mejora la robustez en genes de baja expresión.
Casos de uso:
Preferido para datos con un número moderado a grande de réplicas.
Mejora la detección de genes diferencialmente expresados con bajos conteos.
5.3 pyDESeq2
Método principal:
Es una implementación en Python basada en DESeq2. Internamente utiliza las mismas bases teóricas y estadísticas, ya que está diseñado para ofrecer la misma funcionalidad en un entorno Python.Normalización y estadística:
Reproduce la normalización basada en razones geométricas y los modelos de regresión GLM con dispersión ajustada Bayesiana.
La diferencia clave está en que facilita la integración con bibliotecas y pipelines en Python, como pandas y scikit-learn.
Casos de uso:
Ideal para usuarios que trabajan en ecosistemas Python y no desean depender de R.
Útil en pipelines integrados con otras herramientas de aprendizaje automático o análisis avanzado.
5.4 Resumen
| Aspecto | edgeR | DESeq2 | pyDESeq2 |
|---|---|---|---|
| Base estadística | Binomial negativa + GLM | Binomial negativa + GLM | Igual que DESeq2 |
| Normalización | TMM | Razón geométrica | Igual que DESeq2 |
| Ajuste Bayesiano | Sí (dispersión) | Sí (dispersión y shrinkage) | Igual que DESeq2 |
| Facilidad de uso | R | R | Python |
| Casos recomendados | Pocas réplicas | Réplicas moderadas o altas | Usuarios de Python |
5.5 Resultados obtenidos
Hay que recordar que pyDESeq2 es una librería emergente que trata de implementar la librería original de R DESeq2 en Python. Los resultados muestran una similitud de 86 y 84 genes cuando se comparan las librerías de R DESeq2 y edgeR mientras que se encuentra sólo 1 gen compartido con los resultados de pyDESeq2. Esto mayormente puede deberse a la implementación dada. Sin embargo, no hay que descartar que esten ahí hasta una posterior confirmación (probablemente sólo estén rankeados diferente). Por eso la recomendación es que se sigan usando las librerías de R hasta que se logre hacer una transición completa a Python en términos del análisis transcriptómico.
Reproducibilidad
R
## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.4.2 (2024-10-31 ucrt)
## os Windows 11 x64 (build 26100)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate English_United States.utf8
## ctype English_United States.utf8
## tz America/Mexico_City
## date 2025-01-21
## pandoc 3.2 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-8 2024-09-12 [1] CRAN (R 4.4.1)
## affy 1.84.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## affyio 1.76.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## AnnotationDbi * 1.68.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## apeglm 1.28.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## bbmle 1.0.25.1 2023-12-09 [1] CRAN (R 4.4.2)
## bdsmatrix 1.3-7 2024-03-02 [1] CRAN (R 4.4.0)
## Biobase * 2.66.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## BiocFileCache 2.14.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## BiocGenerics * 0.52.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## BiocIO 1.16.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## BiocManager 1.30.25 2024-08-28 [1] CRAN (R 4.4.2)
## BiocParallel 1.40.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## biomaRt 2.62.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## Biostrings 2.74.1 2024-12-16 [1] Bioconductor 3.20 (R 4.4.2)
## bit 4.5.0.1 2024-12-03 [1] CRAN (R 4.4.2)
## bit64 4.5.2 2024-09-22 [1] CRAN (R 4.4.2)
## bitops 1.0-9 2024-10-03 [1] CRAN (R 4.4.1)
## blob 1.2.4 2023-03-17 [1] CRAN (R 4.4.2)
## cachem 1.1.0 2024-05-16 [1] CRAN (R 4.4.2)
## cli 3.6.3 2024-06-21 [1] CRAN (R 4.4.2)
## coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.4.2)
## codetools 0.2-20 2024-03-31 [2] CRAN (R 4.4.2)
## colorspace 2.1-1 2024-07-26 [1] CRAN (R 4.4.2)
## crayon 1.5.3 2024-06-20 [1] CRAN (R 4.4.2)
## curl 6.0.1 2024-11-14 [1] CRAN (R 4.4.2)
## DBI 1.2.3 2024-06-02 [1] CRAN (R 4.4.2)
## dbplyr 2.5.0 2024-03-19 [1] CRAN (R 4.4.2)
## DelayedArray 0.32.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## DESeq2 * 1.46.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## digest 0.6.37 2024-08-19 [1] CRAN (R 4.4.2)
## dplyr 1.1.4 2023-11-17 [1] CRAN (R 4.4.2)
## edgeR * 4.4.1 2024-12-02 [1] Bioconductor 3.20 (R 4.4.2)
## emdbook 1.3.13 2023-07-03 [1] CRAN (R 4.4.2)
## evaluate 1.0.1 2024-10-10 [1] CRAN (R 4.4.2)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.4.2)
## fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.4.2)
## filelock 1.0.3 2023-12-11 [1] CRAN (R 4.4.2)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.4.2)
## GenomeInfoDb * 1.42.1 2024-11-28 [1] Bioconductor 3.20 (R 4.4.2)
## GenomeInfoDbData 1.2.13 2024-12-22 [1] Bioconductor
## GenomicAlignments 1.42.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## GenomicFeatures * 1.58.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## GenomicRanges * 1.58.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## ggplot2 3.5.1 2024-04-23 [1] CRAN (R 4.4.2)
## ggVennDiagram * 1.5.2 2024-02-20 [1] CRAN (R 4.4.2)
## Glimma * 2.16.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## glue 1.8.0 2024-09-30 [1] CRAN (R 4.4.2)
## GO.db * 3.20.0 2024-12-22 [1] Bioconductor
## graph 1.84.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## gtable 0.3.6 2024-10-25 [1] CRAN (R 4.4.2)
## hexbin 1.28.5 2024-11-13 [1] CRAN (R 4.4.2)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.4.2)
## htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.4.2)
## htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.4.2)
## httr 1.4.7 2023-08-15 [1] CRAN (R 4.4.2)
## httr2 1.0.7 2024-11-26 [1] CRAN (R 4.4.2)
## IRanges * 2.40.1 2024-12-05 [1] Bioconductor 3.20 (R 4.4.2)
## jsonlite 1.8.9 2024-09-20 [1] CRAN (R 4.4.2)
## KEGGREST 1.46.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## KernSmooth 2.23-24 2024-05-17 [2] CRAN (R 4.4.2)
## knitr 1.49 2024-11-08 [1] CRAN (R 4.4.2)
## labeling 0.4.3 2023-08-29 [1] CRAN (R 4.4.0)
## lattice 0.22-6 2024-03-20 [2] CRAN (R 4.4.2)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.4.2)
## limma * 3.62.1 2024-11-03 [1] Bioconductor 3.20 (R 4.4.1)
## locfit 1.5-9.10 2024-06-24 [1] CRAN (R 4.4.2)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.4.2)
## MASS 7.3-61 2024-06-13 [2] CRAN (R 4.4.2)
## Matrix 1.7-1 2024-10-18 [2] CRAN (R 4.4.2)
## MatrixGenerics * 1.18.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## matrixStats * 1.4.1 2024-09-08 [1] CRAN (R 4.4.2)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.4.2)
## munsell 0.5.1 2024-04-01 [1] CRAN (R 4.4.2)
## Mus.musculus * 1.3.1 2024-12-27 [1] Bioconductor
## mvtnorm 1.3-2 2024-11-04 [1] CRAN (R 4.4.2)
## numDeriv 2016.8-1.1 2019-06-06 [1] CRAN (R 4.4.0)
## org.Mm.eg.db * 3.20.0 2024-12-22 [1] Bioconductor
## OrganismDbi * 1.48.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## pheatmap 1.0.12 2019-01-04 [1] CRAN (R 4.4.2)
## pillar 1.10.0 2024-12-17 [1] CRAN (R 4.4.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.4.2)
## plyr 1.8.9 2023-10-02 [1] CRAN (R 4.4.2)
## png 0.1-8 2022-11-29 [1] CRAN (R 4.4.0)
## preprocessCore 1.68.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.4.2)
## progress 1.2.3 2023-12-06 [1] CRAN (R 4.4.2)
## R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.4.0)
## R.oo 1.27.0 2024-11-01 [1] CRAN (R 4.4.1)
## R.utils 2.12.3 2023-11-18 [1] CRAN (R 4.4.2)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.4.2)
## rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.4.2)
## RBGL 1.82.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.4.0)
## Rcpp 1.0.13-1 2024-11-02 [1] CRAN (R 4.4.2)
## RCurl 1.98-1.16 2024-07-11 [1] CRAN (R 4.4.1)
## restfulr 0.0.15 2022-06-16 [1] CRAN (R 4.4.2)
## reticulate * 1.40.0 2024-11-15 [1] CRAN (R 4.4.2)
## rjson 0.2.23 2024-09-16 [1] CRAN (R 4.4.1)
## rlang 1.1.4 2024-06-04 [1] CRAN (R 4.4.2)
## rmarkdown 2.29 2024-11-04 [1] CRAN (R 4.4.2)
## Rsamtools 2.22.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## RSQLite 2.3.9 2024-12-03 [1] CRAN (R 4.4.2)
## rstudioapi 0.17.1 2024-10-22 [1] CRAN (R 4.4.2)
## rtracklayer 1.66.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## S4Arrays 1.6.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## S4Vectors * 0.44.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## scales 1.3.0 2023-11-28 [1] CRAN (R 4.4.2)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.4.2)
## SparseArray 1.6.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## statmod 1.5.0 2023-01-06 [1] CRAN (R 4.4.2)
## stringi 1.8.4 2024-05-06 [1] CRAN (R 4.4.0)
## stringr 1.5.1 2023-11-14 [1] CRAN (R 4.4.2)
## SummarizedExperiment * 1.36.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.4.2)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.4.2)
## TxDb.Mmusculus.UCSC.mm10.knownGene * 3.10.0 2024-12-22 [1] Bioconductor
## txdbmaker 1.2.1 2024-11-25 [1] Bioconductor 3.20 (R 4.4.2)
## UCSC.utils 1.2.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.4.2)
## vsn 3.74.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## withr 3.0.2 2024-10-28 [1] CRAN (R 4.4.2)
## xfun 0.49 2024-10-31 [1] CRAN (R 4.4.2)
## XML 3.99-0.17 2024-06-25 [1] CRAN (R 4.4.1)
## xml2 1.3.6 2023-12-04 [1] CRAN (R 4.4.2)
## XVector 0.46.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
## yaml 2.3.10 2024-07-26 [1] CRAN (R 4.4.1)
## zlibbioc 1.52.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##
## [1] C:/Users/chomb/AppData/Local/R/win-library/4.4
## [2] C:/Program Files/R/R-4.4.2/library
##
## ─ Python configuration ───────────────────────────────────────────────────────
## python: C:/Users/chomb/AppData/Local/Programs/Python/Python312/python.exe
## libpython: C:/Users/chomb/AppData/Local/Programs/Python/Python312/python312.dll
## pythonhome: C:/Users/chomb/AppData/Local/Programs/Python/Python312
## version: 3.12.4 (tags/v3.12.4:8e8a4ba, Jun 6 2024, 19:30:16) [MSC v.1940 64 bit (AMD64)]
## Architecture: 64bit
## numpy: C:/Users/chomb/AppData/Local/Programs/Python/Python312/Lib/site-packages/numpy
## numpy_version: 2.0.2
##
## NOTE: Python version was forced by use_python() function
##
## ──────────────────────────────────────────────────────────────────────────────
Python
## -----
## anndata 0.11.1
## matplotlib 3.10.0
## numpy 2.0.2
## pandas 2.2.3
## pydeseq2 0.4.12
## sanbomics NA
## scanpy 1.10.4
## seaborn 0.13.2
## session_info 1.0.0
## -----
## Python 3.12.4 (tags/v3.12.4:8e8a4ba, Jun 6 2024, 19:30:16) [MSC v.1940 64 bit (AMD64)]
## Windows-11-10.0.26100-SP0
## -----
## Session information updated at 2025-01-21 14:50